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Abstract 



We present a simple, sophisticated method to capture renormalization group 
flow in Monte Carlo simulation, which provides important information of crit- 
ical phenomena. We applied the method to the -D = 3, 4 lattice (j)'^ model and 
obtained a renormalization flow diagram that well reproduces theoretically 
predicted behavior of the continuum (f/^ model. We also show that the method 
can be easily applied to much more complicated models, such as frustrated 
spin models. 

PACS numbers: 02.70.Lq, 75.10.Hk 
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I. INTRODUCTION 



Renormalizat ion-group (RG) theory has drastically improved our perspective about 
phase transition. When it is combined with Monte Carlo (MC) simulations, it becomes a 
very powerful tool to investigate critical phenomena. However, the so-called Monte Carlo 
renormalizat ion group (MCRG) method requires an elaborate scheme and experience, 
and application has been restricted to simple models such as the Ising ferromagnet up 
to now. In this paper we reformulate the MCRG method and present a simple way to 
obtain the RG flow diagram in a MC simulation, which can provide essential information 
about critical phenomena. This paper is organized as follows. In Sec. II, problems in the 
conventional MCRG scheme are explained and remedies for each problem are presented. In 
fact, modification of the conventional scheme leads to use of Binder's parameter [|T^] and the 
second parameter, which is used in recent high-precision numerical analyses ||T5|JT^ . In Sec. 
Ill, expected behavior of RG flow in the D = 4 lattice 0^ model is presented for comparison 
with the MC result. Some caution on MC simulation just at the upper critical dimension 
is also presented. In Sec. IV, details of the MC simulation of the lattice 0^ model are 
described. In Sec. V and VI, results of the MC simulation for D = 3,4 are presented. In 
Sec. VII, we summarize the result and discuss possible applications of the method to more 
complicated models. 



II. MODIFICATION OF THE CONVENTIONAL MCRG METHOD 

In the conventional MCRG scheme critical exponents are calculated from eigen- 

values of the linearized RG transformation: 

dK,{L,b,) 

where Ki{L, b) denotes coupling strength of the z'th term of a block-spin Hamiltonian with 
block size b, and L denotes the original size of the system. 

One problem is that one should use the proper definition of block spin, otherwise Ki{L, b) 
goes to zero or infinity as b becomes large even at the critical point. For the Ising model, 
the majority rule is one of the options. However, it cannot be extended to more complicated 
models. There are infinite kinds of definitions and location of the fixed point (along an 
irrelevant direction) depends on it 0. 

In this paper, we use one that seems most simple and suitable for Monte Carlo simula- 
tion, 

Sb = Sb/^< si >, (2) 

where Sb denotes summation of all spins in the block. This definition is implicitly used in the 
definition of Binder's parameter. Actually, we use coarse-grained spin obtained by cutting 
high-momentum modes off, as in Ref . . However, we set 6 = L at the last stage, therefore 
real- space RG and momentum-space RG do not make much difference. 

Another problem is that the block size b should be much smaller than L, otherwise the 
behavior of the block spin is affected by a boundary condition that is nonuniversal. Thus one 
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faces a tradeoff that for larger b, behavior of block spins become more asymptotic (couplings 
of irrelevant terms become smaller) but their behavior deviates from bulk one. The solution 
is to use the following matrix 



dKi{bL,bL) _dKi{bL,bL) dK,{bL,b) 
dKj{L,L) ~ dKj{bL,b) ' dKj{bL,l) 



dKi{L,L) 
dK,{L,l)\ 



1 -1 



(3) 



instead of (|T|). At the critical point, 

dKi{bL,bL) 



dKiiL,L) 



dKj{bL,b) dKj{L,l) 



(4) 



should be satisfied and eigenvalues of the matrix (|^) coincide with that of (|l]). Thus one can 
use as large a block size as the system size. 

Yet another, and most severe problem, is referred to as the "redundancy problem" 
In the conventional MCRG scheme, one should observe very large numbers of block-spin 
interaction terms to construct a fixed-point Hamiltonian. However, in the continuum (j)"^ 
theory, there are only two kinds of independent coupling constants (mass and coupling). 
The lattice Hamiltonian approaches this continuum asymptotically after the momentum- 
space RG transformation, and the two terms are enough to observe essential RG flow. 

Let us consider the Hamiltonian deflned on D-dimensional continuum space: 



n 



I [V0(x)]^ + |c 



(x)2 + /50(x) 



(5) 



Note that it is not a microscopic Hamiltonian: it is a phenomenological Hamiltonian and 
the parameters 7, a, (3 should be determined to reproduce experimental results. In other 
words, it is the logarithm of probability of observing a speciflc conflguration of some physical 
quantity 0(x) in the experiment. This means that short length scale fluctuations beneath 
the resolution / of the observation device are already integrated out and absorbed into the 
parameters 7,a,/3. Thus parameters depend on the cutoff length / and should be denoted 
by jijaijjSi. We also denote the physical quantity 0(x) averaged over a volume of linear 
length / by 0i(x). 

We replace the parameters in (W) by "regularized" ones deflned as follows: 



ai = l^ < (t)i{xf > ai, Pi = l'' < 0/(x)' pi, 7z = Z^"' < 
In terms of this regularization, the Hamiltonian is expressed as 



jD-2 



'|pi|<7r 

where 0(p) = 0(p/) < (pf 



a -\- 7p^ - - 
rfp 0(p)0(-p) 

,p|<7r ^ 
c/pirfp2(ip3(ip45(XlPi)0(Pl)0(P2)0(P3)0(P4), 



L 



(6) 



> 



-1/2 



All possible values of the regularized parameter fall on 



a two-dimensional manifold on which < (j)i{x) >= 1 is satisfled. This regularization is 
suitable for MC simulations, compared to the fleld theoretical one, 7 = 1. Figure ^ and 
Fig.p| show schematic RG flow of these regularized parameters in D = 3 and 4, respectively: 
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high temperature, low temperature, Gaussian, and Wilson-Fisher fixed point are denoted 
by H,L,G, and WF, respectively. 

Now let us return to the lattice model defined as follows: 

^L = |E(^^-^.r + Ef^? + /54- (7) 

<ij> i 

where the summation J2<ij> suns over all nearest neighbor pairs. If we apply the momentum- 
space RG transform of factor 6 to (^, the renormalized Hamiltonian takes the following 
form: 

-(T u\ (L/b) 

H{L, b) = ^^{L/b)-^ s-,{p)s,{-p) (8) 

|p|<i" 

-(T U\ (L/b) 

+ ^^{L/by^ E p^.-,(p).-,(-p) (9) 

|P|<7I' 

(L/b) 

+ ^{L,b){L/b)-^'' E s,(pi)s,(p2)s,(p3)s,(p4), (10) 

Pl+P2+P3+P4=0 

where Sb{p) = s{bp) < sl{x) >~^/2 and the symbol I]p^^ denotes summation over 
0, ±27r/A^, ±47r/A^, ■ ■ ■ , for each component of p, s{p) denotes Fourier components of Sj, 
and 

.,(x) = Y: e^^-.(p) (11) 

|p|<7r/b 

is coarse-grained spin at x. Actually, higher terms such as 0{sb^) are present in H{L, b) but 
they vanish as L and b become large. 

Note that < X]p ■5b(p)sf,(— p) >= const by definition and there are only two kinds of 
interesting terms. If we set b = L, the term (p!oD is exactly the Binder's parameter and we 
denote it by Bl- The term (H) vanishes if we set b = L, so we must stop RG transform at 
b = L/2. Then it becomes 

27r^ < s(kOs(-ki) > 
2-^ < s(0)2 + 2s(ki)s(-ki) >' ^ 

where ki = (27r/L, 0, ■ ■ ■ , 0). For simplicity, we use the following quantity: 

_ < g(ki)s(-ki) > 

= < sm > ■ ■ ^ ^ 

Recently, several different parameters have been proposed as the second parameter Cl to 
estimate and eliminate the sub leading scaling field. We will review these works in Sec. II A. 
Now consider the following matrix: 



d{BbL,CbL) _ d{BbL,CbL) ( d{BL,CL) 
d{BL,CL) 9(A,7i) 'U(A,7i) 



(14) 



instead of 



d[^{bL,bL),^{bL,bL/2)] 



At the fixed point, 



d[PiL,L),jiL,L/2)] 



(15) 



(16) 



(17) 



is satisfied and eigenvalues of d{BhL, Cijl) / di^Bi, Cl) coincide with that of d{(3i,, %)/d{f3i, 71) 
as long as 9(5^, C£)/9(/5i, 71) is a nonsingular matrix. 

Thus scaling behavior of renormalized parameters can be extracted from that of {Bl, Cl), 
as long as \d{BL,CL)/d{/3i,'yi)\ 7^ 0: If we draw arrows fYom{BL,CL) to {BbL,CbL) in the 
Bl-Cl plane, the renormalization flow diagram of factor b is obtained. From this diagram, 
one can determine whether the MC result is asymptotic enough or not, by checking whether 
{Bl, Cl) converge to a fixed point. When a subleading scaling field is expected to be very 



large, such as in the recent Monte Carlo studies of five-dimensional (5D) Ising model pA|-p!3 
this RG flow diagram is very useful compared to a many parameter fit to a single observable 
Bl- 

A linearized RGT matrix = d(BbL, ChL) /d{BL, Cl) is calculated, for example, from 
linear fitting 

( BbL{a,(3,-f) \ _-^f BL{a,(3,-f) \ / Bo\ ^^g^ 



\C,L{a,(3,^) ) ^'\CL{a,(3,^) 

where Rf, and {Bq,Co) are fitting parameters and we use values of {Bl,Cl) and {BbL,CbL) 
at several different parameters 7 near the fixed point. Selection of this parameter is a 
delicate problem: when it is too close to the fixed point, {Bl, Cl) and {BhL, C^l) become 
very close to each other and the RG flow is buried in statistical errors, while when it is 
far from the fixed point, nonlinear dependence of {Bl, Cl) on {a, {3, 7) induces a systematic 
error. Thus the parameter range for the fitting should be determined carefully. 

Of course Rj, can be calculated from a, j3, 7 derivatives of Bl and Cl at the fixed point. 
However, derivatives with respect to irrelevant direction vanish as L becomes large and 
buried in statistical errors. Thus values of Bl, Cl at parameters well apart from the fixed 
point along an irrelevant direction are needed to calculate the second eigenvalue. 



A. The second parameter in literature 

Here we compare our definition of Cl, "the second parameter" (the first one is the 
Binder's parameter), with preceding works. 

Ballesteros et al. [|14| used a finite-size correlation length defined below as 

_ < > I < 0(ki)0(-ki) > -1 

^ - sin^dkil) ' ^'^^ 
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which is related to Cl as sin^(|ki|) = — 1 + 1/Cl- 



Hasenbusch |]15[ used the ratio between the partition function of a periodic and an- 
tiperiodic Hamiltonian. In momentum representation, the Hamiltonian for an antiperiodic 
boundary condition (APBC) is obtained by replacing k summation of the lattice Hamilto- 
nian by fc^ = ±7r/L, ±37r/L, ■ ■ ■ for each direction to which the APBC is imposed. Let 
us denote the partition function of a periodic and antiperiodic Hamiltonian by Zp and Za, 
respectively. It then reads: 

ZJZ, = {IT^Yl'l (20) 

_ /P0exp(-gp) eMHp - Hg) _ ^ 

TTfTi 7 — Tr\ ~< exp(iip - Ha) >p [21) 

where Hp and Ha denote the periodic and antiperiodic Hamiltonian, respectively, < ■ ■ ■ >p 
denotes average with respect to Hp. When the APBC is imposed to a direction of ei, it 
reads 

Ha-Hpr^Y. ^0(k)0(-k) (22) 

for large L. One can see that Za/Zp has a similar form as that of Cl- In the real-space 
RG scheme, both Cl and Za/Zp can be regarded as an effective coupling between two block 
spins defined on L x L x ■ ■ ■ x L/2 block. 



III. PERTURBATION EXPANSION AT D = A 

Near the Gaussian fixed point, finite-size behavior of Bl and Cl can be predicted 



from finite-size perturbation theory proposed by Chen and Dohm [jT6[. Note that, when 
D = 4, there is one kind of divergent subdiagram (Fig. ^ whose factor is proportional 
to [/?i(Co + Cl lnL)/7^] at the critical region, where Co,Ci > are some constants. Thus 
perturbation is restricted to the range /?i(Co + Ci lnL)/7^ <^ 1 and one cannot set L — > cxd. 
However, the limit for L rapidly diverges as we approach the Gaussian fixed point, and good 
agreement between perturbation theory and Monte Carlo data is expected for a certain pa- 
rameter range and lattice size. Then one can predict scaling behavior for large L, which 
is far beyond the computational limit of Monte Carlo simulation, from the perturbation 
theory. 

Here we investigate finite-size behavior of Bl and Cl at the finite-size critical point (to 
one-loop order): 

where J(k) = J2fi{^ — cosk^). 

Bl takes the so-called zero-mode value: 
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^™ - . :r'jj :r - 2.1844. (24) 



/ d% exp(-$^)$^/ / d% exp(-$ ^ 
[/ d^o exp(-<|.4)$2/ / d$o exp(-$^)] 

As for Cl, to one- loop order, 



For large L, X]k7^o[2^(k)]~^ ~ + 742lnL where Ai and are some positive constant. 
Thus the plot of {Bl,Cl) for a certain range of parameters including critical temperature 
approaches (-Bzm, 0) as L increases or /? decreases, indicating that the Gaussian fixed point 
is infrared stable. However, approach to the point {Bzm,0) for increasing L is extremely 
slow, and one cannot expect asymptotic Gaussian behavior in MC simulations, even when 
L is very large. Critical exponents estimated from MC results may also differ from the 
Gaussian (classical) value and depend on the bare parameters. Thus one can extract only 
very restricted information from MC simulations at the upper critical dimension. 



IV. DETAIL OF MONTE CARLO SIMULATION 

We investigated the L x L x L system with L = 8, 16, 32 for D = 3 and the L x L x 
L X L system with L = 4, 8, 16 for D = 4, imposing a periodic boundary condition on each 
directions. 

We used the following Hamiltonian: 

^ = ^E(^^-^.r + fE^? + /3E4 (26) 
^ <ij> ^ i i 

where —oc < Si < oo denotes the spin on the site i. Since we cannot use the regularization 
condition < >= 1 before the simulation, we used the following one 

/^0exp(-f0^-/j0^)0^ _ 
/rf0exp(-f02_;504) 

We used several fixed a, (3 and tuned 7 to reach the critical region. Actual values of param- 
eters are listed in the later sections. 

For parameters well apart from the Gaussian fixed point, 4L^~^ single cluster flips |jT7 



and 16 Metropolis sweeps are performed between successive observations. In the cluster- 
update stage, length of Sj is kept fixed and only its sign is changed. In the Metropolis update 
step, a new value for Si is chosen uniformly from a range exp(— as^ — /9s^) > exp(— 3.0). By 
the cluster update, all spins are flipped four times on average between observations. For all 
observed quantities, the correlation coefficient between successively observed values was less 
than 0.2. 

As we approach the Gaussian fixed point, parameters behave as follows: 

a ~ 0, /? ~ const, 7 — >■ 00. (27) 

Thus the nearest neighbor coupling 'jSiSj/2 tends to diverge. Since the bond-cutting prob- 
ability of Wolff's algorithm is exp(— 7SjSj), the cluster update tends to end up with flipping 
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the whole system and does not accelerate the simulation. In these cases we increased the 
number of Metropolis sweeps until the above mentioned condition is satisfied. 

Thermal averages of observables at 7 slightly away from the actually used value were 
calculated using reweighting techniques |T^. For all system sizes, at least 0.8 x 10^ observa- 
tions were done after thermalization, and statistical errors were estimated by the jack-knife 
procedure. Multiplicative lagged Fibonacci sequence Rt = Rt-mss x Rt-4i87 (mod2^^) was 
used as a random number generator. All runs were performed on VPP-300 at JAERI. 



V. RESULT FOR D = 3 

Simulations were performed at the following parameters (see also Fig. P for L = 8, 16, 32; 



Ising case: (j) 


(x) = ±1, 


7 = 


0.2217; 






Case A: 


a/2 


= -2.6159, 


[3 


= 1.0948, 


7 


= 0.3040; 


Case B: 


a/2 


= -1.6655, 




= 0.6935, 


7 


= 0.3545; 


Case C: 


a/2 


= -0.8786, 




= 0.3938, 


7 


= 0.4700; 


Case D: 


a/2 


= -0.5209, 


/3 


= 0.2713, 


7 


= 0.6260 



Figure § shows the RG fiow diagram of Bl and Cl near the Wilson-Fisher RG fixed 
point. All lines are drawn from {Bl,Cl) to {B2l,C2l)- Dashed and solid lines correspond 
to L = 8 and L = 16, respectively. One can see that the RG fiow shown in Fig. |I] is well 
reproduced. Moreover, the linear fitting procedure ([18|) provides an estimate for critical 
exponents as z/ = 0.69(3), tu = 0.74(10) from L = 8, 16 data and u = 0.653(10), u; = 0.7(2) 
from L = 16, 32 data, which is in agreement with the most recent e-expansion result u = 
0.6305(25), = 0.814(10) ^ and the Monte Carlo result u = 0.6296(3), 0^ = 0.845(10) |T§. 
Thus one can see that the two observables Bl and Cl are enough to capture the essential 
RG fiow. 



VI. RESULT FOR D = 4 

Simulations were performed at the following parameters, 
Ising case: 0(x) = ±1, 71 = 0.1495; 
Case A: a/2 = -0.3226, p = 0.2082, 7 = 0.54; 
Case B: a/2 = -0.1383, /5 = 0.1530, 7 = 1.0 

for L = 4, 8, and 16. Near the Gaussian fixed point (Case B), there is a severe, critical 
slowing down (owing to large fiuctuation of the k = mode), which cannot be removed by 
the cluster update, and we did not perform L = 16 simulation. For Case B, perturbation 
expansion agrees well with the MC result: Fig.^a) and H(b) show the plot of Bl and Cl, 
respectively, against 71 for fixed (a,/3), together with perturbation results. Note that there 
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are no free parameters to be fitted, unlike tlie Ising case and the agreement is both 
quahtative and quantitative. 

Figure |^ shows a RG flow diagram of and Cl obtained from MC simulations. All 
lines are drawn from {Bl,Cl) to {B2l,C2l)- Dashed and solid lines correspond to L = 4 
and L = 8, respectively. Simulations at a parameter closer to the Gaussian fixed point than 
Case B are very difficult owing to aforementioned critical slowing down. Instead, finite-size 
perturbation provides reliable results near the Gaussian fixed point and it indicates that the 
plot of {Bl, Cl) approaches the Gaussian fixed point as L increases. Thus one can conclude 
that there is no RG fixed point except for the infrared-stable Gaussian fixed point. 



VII. CONCLUSION 

The renormalization-group flow diagram obtained by the method presented in this paper 
provides qualitative information such as the stability of a specific RG fixed point against 
some perturbation, as well as quantitative improvement of the estimated value of a critical 
exponent by eliminating leading correction-to-scaling terms. However, in lattice models, 
there exists 0{1/ L) systematic error owing to the substitution of integral by finite summation 
oil/L mesh, and one cannot get rid of this as long as the finite lattice system is concerned. 

Our method can be easily extended to ip^ models with several or unusual quartic coupling 
constant (s): 

H = I dx{V(t)f + + ^ I3n E Cijki{n)(t)i(t)jMu (28) 

n ijkl 



such as the chiral 0{2n) model of a triangular antiferromagnet and the Ginzburg 



Landau model of a type-II superconductor under a weak or strong magnetic field, with or 
without point / columnar impurities (see p2| , p3| for transition of a pure system under a strong 
field). The multiple quartic term tends to generate an irrelevant operator whose correction 
exponent is very small, making it difficult to observe asymptotic behavior in MC simulations. 
Thus critical behavior of these models has been a controversial issue and application of the 
MCRG method to these models seems very interesting. A RG flow diagram of regularized 
quartic terms < Cijki{n)(f)i(j)j(j)k4>i > / < J2i (pf will reveal the critical behavior, as in Ref. 
of the Q = 4 antiferromagnetic Potts model. 

Another important problem is the estimation of the correction exponent for the 0^ term 
at the WF fixed point. An e expansion analysis indicates that it becomes positive at the 
WF fixed point [Q. However, whether the exponent is larger or smaller than u ^ 0.8 
of (j)^ theory, which we assumed as a leading correction, should be confirmed numerically. 
Similarly, the effect of sixfold anisotropy on the critical behavior of a three-dimensional XY 
model is another interesting subject since the anisotropy is expressed by the 0(0^) term. 
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FIGURES 
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FIG. 1. Renormalization flow of regularized parameter in £) = 3. 




FIG. 2. Renormalization flow of regularized parameter in = 4. 




FIG. 3. Divergent subdiagram in D = 4. 
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FIG. 4. Parameters at which simulations were performed. 
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FIG. 5. RG flow near the Wilson-Fisher fixed point in L> = 3. All lines are drawn from {Bl, Cl) 

to {B2L,C2l). 
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FIG. 6. Plot of (a) Bl and (b) Cl against 71 for the Case B in D = 4. "MC" denotes Monte 
Carlo result and "theory" denotes finite-size perturbation result. 
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FIG. 7. RG flow of Bl and Cl , obtained by MC simulation in D = 4. All lines are drawn 
from {Bl,Cl) to {B2l,C2l)- Dashed and solid lines correspond to L = 4 and L = 8, respectively. 



13 



